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We present an equation-free dynamic renormalization approach to the computational study of 
coarse-grained, self-similar dynamic behavior in multidimensional particle systems. The approach is 
aimed at problems for which evolution equations for coarse-scale observables (e.g. particle density) 
are not explicitly available. Our illustrative example involves Brownian particles in a 2D Couette 
flow; marginal and conditional Inverse Cumulative Distribution Functions (ICDFs) constitute the 
macroscopic observables of the evolving particle distributions. 
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I. INTRODUCTION 

Multiscalc phenomena arise naturally in science and 
engineering. In many cases of current research inter- 
est, physical models are available at a fine, microscopic 
scale (atomistic, stochastic, agent-based), while we want 
to study the system behavior at a coarse-grained, macro- 
scopic level. Macroscopic, coarse-grained evolution equa- 
tions, even when they conceptually exist are often un- 
available in closed form, due to the lack of accurate ex- 
plicit closures. The Equation-Free computational frame- 
work ([1,0, 0,0) nas been recently proposed for the 
computer-assisted study of precisely such complex, mul- 
tiscale problems: those for whose macroscopic behavior 
no explicit coarse-grained evolution equations are avail- 
able. Equation-free methods utilize so-called coarse time- 
steppers, which are used to numerically analyze coarse- 
grained behavior through appropriately designed short 
computational experiments performed by the fine-scale 
models; coarse timestepping is closely related to optimal 
prediction in the work of Chorin and coworkers [|| Q 
(see also the discussion in 0] ) . Quantities necessary for 
traditional continuum numerical analysis (residuals, the 
action of Jacobians) are estimated on demand from these 
short fine-scale computational experiments, rather than 
evaluated from closed-form macrosopic equations. This 
computational "enabling technology" has been used to 
perform integration, fixed point computation, numerical 
stability and bifurcation analysis as well as control and 
optimization for the coarse-grained behavior of a num- 
ber of fine-scale model types (Lattice Boltzmann, Monte 
Carlo, Molecular Dynamics, Brownian Dynamics etc., see 
the references in 0, 0] ) . 

For problems whose macroscopic behavior is character- 
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ized by scale invariance, dynamic renormalization meth- 
ods 0, E3 provide tools for locating self-similar solu- 
tions and their scaling exponents 0, 13 • Recently, we 
have combined equation-free computation with dynamic 
renormalization to obtain coarse-grained self-similar so- 
lutions based on direct microscopic simulation 0, Q; 
our particular implementation used a template-based ap- 
proach 0, HHli. Working in a frame of refer- 
ence that expands (or shrinks) along with the macro- 
scopic system observables turns the self-similar problem 
into a steady state one, for which fixed point techniques 
can be applied. Equation free dynamic renormalization 
techniques are used here to investigate macroscopic self- 
similarity in multidimensional particle system dynam- 
ics. This is made possible through a coarse time-stepper 
which utilizes marginal and one-dimensional conditional 
inverse cumulative distribution functions (ICDFs) as the 
coarse observables of multidimensional particle systems. 
Our illustrative example consists of the (self-similar) dy- 
namics of Brownian particles in a Couette flow; the re- 
sults are validated using known analytical solutions. The 
paper is organized as follows: we start with a brief de- 
scription of the coarse time-stepper in our framework. 
We then discuss equation-free dynamic renormalization, 
present our illustrative example, demonstrate the compu- 
tation of its self-similar shapes and exponents, and con- 
clude with a brief discussion of limitations and possible 
extensions of the approach. 



II. A COARSE TIME-STEPPER FOR 
MULTIDIMENSIONAL PARTICLE SYSTEMS 

In the equation- free framework, short, appropriately 
initialized computational experiments with a fine scale 
model are used to construct the coarse time-stepper- 
the basic element capturing the dynamical interaction 
between coarse-scale observables and fine-scale model 



2 



Coarse state 1 



Coarse state 2 



Lifting 



Restriction 



Fine state 1 



Fine-scale evolution 



Fine state 2 



FIG. 1: Schematic of the coarse time-stepper 



states. There are essentially three components in a coarse 
time-stepper, namely, lifting, fine-scale evolution and re- 
striction |8| (Fig. Lifting is a transformation that con- 
verts coarse-scale observables to one or more consistent 
fine-scale realizations; restriction is the transformation 
in the reverse direction, from fine scale states to coarse 
observables. The specific manner in which these lifting 
and restriction operators are implemented yields differ- 
ent inter-scale exchange protocols; one must test that the 
results at the macroscopic level are relatively insensitive 
to small details in the protocols (see the discussion in 
0). 

In macroscopically multidimensional particle systems, 
particle positions are a natural component of the fine 
scale state, and marginal and conditional ICDFs of parti- 
cle positions are good candidate coarse-scale observables. 
Assuming smoothness, a finite number of marginal and 
conditional ICDFs can be used to recover the distribu- 
tion of particle positions (e.g. through interpolation). 
Our coarse observables are the leading coefficients of the 
projection of these ICDFs on an appropriate orthonormal 
(Shifted Legendre) polynomial basis. Particle positions in 
one direction are generated using the marginal ICDF, and 
the corresponding particle positions in the second direc- 
tion are generated using the conditional ICDFs in our 
lifting step. For the restriction step, after constructing a 
smooth multidimensional CDF from particle positions by 
simple interpolation, marginal and conditional CDFs are 
first obtained, and then their inverse CDFs are interpo- 
lated and their leading expansion coefficients computed. 



III. EQUATION-FREE DYNAMIC 
RENORMALIZATION 



For many macroscale systems of practical interest, if 
the PDEs describing their evolution are scale invariant, 
they may possess self-similar solutions 0, ITlj . Dynamic 
renormalization procedures have been used for the study 
of such self-similar solutions [H El El; recently the 
template-based method for studying the dynamics of 
problems with symmetry [22^ has been extended to study 
the dynamics of problems with scale invariance |lL l2lLl23| . 
When scale invariant evolution PDEs are explicitly avail- 
able, a template based approach can be used to derive 
dynamical equations (termed "MN-dynamics" ) for the 
renormalized self-similar solutions and similarity expo- 



nents EI- The idea underlying this approach, especially 
that of employing template conditions, can be used to 
obtain renormalized self-similar solutions and similarity 
exponents in systems whose macroscale governing equa- 
tions are not explicitly available Q (see also the approach 



m0). 

Consider a PDE of the form 



dF 
~dt 



(1) 



where F(x,y,i) is a CDF of particle positions. We as- 
sume that for the differential operator D xy there exist 
constants p and a such that 



D xy (f(j,^))=A a D uv (f(u,v)), u 



£ , _ v 

A' V ~ Ap 
(2) 

for any real function /, real value A > and coordi- 
nate (a;, y) (there is no amplitude rcscaling since this is 
a CDF). If a self-similar solution F(x,y,t) exists, it can 
be written as 



F(x,y,t) = U( 



tJ 



{cs) a ' {cs) a P 



;c), 



(3) 



where c is a constant parametrizing the family of renor- 
malized shapes, s = t — to, (to is the blowup time for 
problems with finite time singularities, whether forward 
or backward in time). 
Therefore, 

aa = — 1, (4) 

and U satisfies the PDE, 

-auU u - apvU v = c~ 1 D uv (U), (5) 

where u = x/(cs) a , v = y/(cs) ap . For D xy satisfying 
Eqn.©, the constant a is determined by D xy itself and 
the similarity exponent a can be obtained by Eqn.QJ. 

If the equation is not explicitly available, one cannot 
analytically obtain the exponents a and p; tests have to 
be devised for finding these constants -and thus testing 
the scale invariance of the operator- before we embark 
upon the computation of the self-similar solutions them- 
selves. 

For D xy satisfying the (unknown) equation ©, the 
constants p and a can be obtained using a black box 
simulator of the equation as follows: Since the unknown 
Eqn. J3J) is valid for any real function /, let / be a test 
function, and let us choose the two points {u\,v\) and 
(u 2 ,v 2 ) and let (xi,yi) = (u\A,viA p ) and (2:2,1/2) = 
(112A, V2A P ), where A is arbitrarily chosen as a positive 
real value. Then the following two relations should hold 
for essentially arbitrary test functions f(x,y) and points 
(ui,Vi),i = 1,2: 

D xy (f(j,j^))(x u yi) = A a D uv (f(u,v))(u uVl ), 
D Mj^))(^,m) = A a D uv (f(u,v))(u 2 ,v 2 ).(6) 
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Comparing the above two equations, we have 
D xy(f(% ,%;)){xi,yi) _ D uv (f(u,v))( Ul ,vi) 



D xy{f{% $?)){x2,y2) D uv (f(u, V))(U 2 , V 2 ) 



(7) 



Eqn.Q is solved (using Newton's method) for the 
constant p. When D xy is not explicitly available, 
D xy (f(x,y)) can be obtained by running the micro- 
simulator for a short time and numerically estimating 
the derivative Given p, the constant a is calculated 
by 



A D uv (f(u,v))(ui,v 1 ) 



(8) 



Clearly, other test functions, as well as conditions eval- 
uated at other points (or nonlocal versions of the condi- 
tions using templates) can be used; care must be taken 
also to ensure the finiteness of the estimated quantities. 

Given a and p, to determine the self-similar shape of 
the solution, we consider the general scaling 



F{x,y,t) = lu(- 



y 



(9) 



where A(t) is an unknown function. The PDE becomes 

L0 t - ^uuj u - ^^vuj v = A a D uv (uj). (10) 

Evidently, U and u> are both rescaled CDFs. 

A single template condition coupled with equation itTUl) 
is needed to solve for both u){u,v,t) and A{t) at every 
time step. In this paper we choose this template condi- 
tion to be 

w(e, oo, t) = m, e < 0, < to < 0.5, (11) 

where e and to are both (essentially arbitrary) constants. 
This template condition has the following physical mean- 
ing: the u-coordinates corresponding to the rescaled 
marginal CDF ujjj = m have the same value e for all 
t. 

Applying this template to Ean. H10|) and assuming 
jj£(e,v,t) decays exponentially as v — > oo, we have 



A t e^-{e,oo,t) + A a+1 D uv (u)(e,oo,t) = 0. 
ou 



(12) 



Equations (|10|) and i|12|) are solved for the rescaled 
CDF oj(u, v, t) and rescaling variable A(t) if the operator 
D xy is explicitly known. As the time t — > oo, tu may ap- 
proach a steady state, which is a stable self-similar shape 
for the solutions to Eqn.JTJ. 

The value for a can be calculated in the long-time limit 
(i.e., after u> reaches steady state). Indeed, let t\ and t% 
be distinct times after w reaches the steady state, then 
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FIG. 2: Illustration for the coarse renormalization 

In many cases the macroscale equation for the CDF of 
particle positions may not be explicitly known. However, 
the template conditions can still be used to rcnormalizc 
the CDF itself, evolved via microscale models, and rescal- 
ing variables are obtained during the course of renormal- 
ization. In these cases, a coarse time-stepper is employed 
to evolve the coarse system using the the marginal and 
conditional ICDFs (actually, the leading coefficients of 
their projection on an appropriate basis) as coarse ob- 
scrvables. The procedure for equation-free renormaliza- 
tion, depicted in Fig. consists of the following six 
steps: 

1. Generate the marginal and conditional ICDFs ac- 
cording to the initial CDF or according to the co- 
efficients of dominant modes of the ICDFs. 

2. Generate particle positions based on these ICDFs 
using the lifting procedure in the coarse time- 
stepper. 

3. Evolve particle positions over a time interval T' us- 
ing the fine-scale model. 

4. Obtain ICDFs from particle positions using the re- 
striction procedure in the coarse time-stepper. 

5. Rescale the marginal ICDF according to the tem- 
plate condition and obtain the rescaling variable 
A. We then rescale the conditional ICDFs by A p 
(recall that p has been independently computed 
through Eqn. J7J). This step can be justified 
by Eqn.@. Indeed, obtaining the rescaled solu- 
tion LUk+i from cuk via the dynamics (fTUf) and jTJjl 
is equivalent to starting from the initial condition 
ujk, running the original dynamics for a while to 
get Ffe+i, obtaining the rescaling variable A, and 
rescaling Fk+i by A and A p respectively in the x 
and y directions. 

6. Project the rescaled ICDFs onto the appropriate 
orthonormal basis and obtain the coefficients of its 
dominant modes. Go back to step 1. 

The above procedure can be viewed as an iterative al- 
gorithm to solve for the fixed point of a nonlinear oper- 
ator $t', written as, 



A t (t 2 ) A t {t t ) 



a - $t'(q) = . 



(14) 
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This fixed point can be written in component form as, 
a r i>p ,i = l,--- ,M + l,p = 0,-.- ,Pi,or 

a = (a^ ; a i,i' ' ' ' ' a i,P' a 2,o> a 2,i' • ' ' ! a 2,p ; 

' ' ' ' a M+l,0i a M+l,l! " " ) a M+l,p) ! 

where •, j = 0, 1, ■ ■ • , P stands for the projection of 
the marginal ICDF onto a j'th-order mode, a^-,* = 
2, ■ • • , M+ 1; j — 0, 1, ■ • • , P the projection of the % — 1th 
conditional ICDF onto a jth-order mode, and the super- 
script r refers to the fact that this is in fact the renor- 
malized self-similar shape. These coefficients correspond 
to the renormalized ICDFs and CDF of the multidimen- 
sional particle system. Equation l|14|) can also be solved 
using any numerical algorithm such as direct iteration 
or matrix-free implementations of Newton's method (e.g. 
Newton-GMRES, H3). 



In the above equation, the operator D xy is written as 

d f x d D 2 d 2 
D ^ = - X 8-y + Lo-y dX + ^- (19) 

This operator satisfies the invariance property 
with constants p = 3 and a = — 2. 

The analytical self-similar solution to the PDF equa- 
tion JT3 is 

/q _( 6(y-0.5x(t-t )) 2 2 \ 

P XY (x,y,t) = -v* e~\ D ^ + ^^) , 
TrD 2 (t - t ) 2 

(20) 

where to is the blowup time (forward or backward in 
time); the corresponding CDF self-similar solution to 
(HHJ is 

Fxy{x,y,t) = 



IV. BROWNIAN PARTICLES IN COUETTE 
FLOW 

We now illustrate our computational approach to self- 
similarity in a two-dimensional Brownian model of par- 
ticle dispersion in Couette flow [2(|. Let X(t) and Y(t) 
represent particle positions in the x and y directions re- 
spectively at time t on the plane. The particle positions 
evolve according to 



dX(t) = Ddu> x (t), dY(t) = Xdt, 



(15) 



where u>x if) is a Wiener process Q and D is the diffusion 
coefficient. The discretized dynamics of (|15|l is given by 

m 

X k+1 =X k + Drj X ,kVM, Y k+ i = Y k + X k dt, (16) 

where rjx,k are i.i.d. standard Gaussian random vari- 
ables. 

The dynamics l|15(l represent the motion of particles 
which only diffuse in the x-direction in a Couette flow 
(see the Discussion for both x and y diffusion). It can 
be shown that the coarse-scale dynamics for the PDF, 
Pxy(x, y, t), of a particle position, corresponding to the 
fine-scale dynamics <|15|) . is governed by the following 
equation 



dy 

m 



dPxr dPxY D 2 d 2 P XY 



dt 



dy 



dx 2 



(17) 



where Pxy is assumed to be 2nd-ordcr diffcrcntiablc. 
Hence the dynamics for the CDF, Fxy{x, y, t), associ- 
ated with (|T7jl is given by 



dF 



XY 



dF 



XY 



dt 



dy 



dF 



XY 



dy 



dx 



D 2 d 2 F XY 
2 dx 2 



(18) 



V3 



irD 2 {t-t ) 2 J_ 



,oiy-u.oatt«-« ;; i x" \ 

e o^'-'o) 3 wit-tor dydz. 



oo J — oo 



Let u' = (e(f _ t a o)) i/2 and v> = (c(t _f o))3/2 ; then 
F XY {x,y,t) =F uv (u',v') = 



(21) 



V3 



irD 2 /c 2 



- 



'wT^'dvdu. (22) 



Hence for the intcgro-diffcrcntial equation l|18|) , the sim- 
ilarity exponent a in |J3J| is a = 1/2. For the CDF 
in 122|) . its std.'s in two directions and correlation arc 
a x = D/c 1 / 2 , a Y = D/(V3c 3 / 2 ), and p XY = V3/2, re- 
spectively. 



V. NUMERICAL RESULTS 

We now compute the self-similar solutions without the 
macroscopic equations, based only on the particle sim- 
ulator. The invariance property of the (presumed un- 
available) macroscale differential operator D xy has to be 
established first. Then the fixed-point algorithm is used 
to solve for the renormalized CDF steady-state shape, 
and the similarity exponent is computed. The diffusion 
coefficient D and simulation time step At in the fine-scale 
model are set to S.Ocm/s 1 / 2 and 0.01s, respectively. An 
ensemble of 9000 (N = 9000) particles is typically used 
in the fine-scale simulations. 

Without the macroscale equation JTSJ, we use New- 
ton's method to solve equation for p. Our test 
function / is a 2-dimensional joint Gaussian, f(x,y) = 
l/4.5 2 A r (a;/4.5)A r (?;/4.5), where iV(a;) andiV(y) are stan- 
dard Gaussian distributions. The positive real number 
A = 2.0. The two coordinates (ui,Vi) and (1*2,^2) in 
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No. of Iterations 


P 


a 





6.0 


-3.86287 


1 


2.79441 


-1.83047 


2 


2.97849 


-1.97927 


3 


3.00132 


-2.00114 


4 


2.99804 


-1.99732 


5 


2.99670 


-1.99612 


6 


2.99441 


-1.99675 


7 


2.99343 


-1.99113 


8 


2.99400 


-1.99459 



TABLE I: Iterations in the computation of the constants p 
and a 



are chosen as (—2.5, —2.5) and (3.5, 3.5), respectively. To 
reduce fluctuations of values for the operator D xy , 200 
copies of values for D xy and D uv are averaged in the 
computation. 

Starting from the initial value p = 6.0, iterative values 
for p are stabilized at 3.0 after 3 iterations. Accordingly, 
the converged value for a is approximately —2.0 (TableP) 
(other test functions and test point selections gave similar 
results) . 

This suggests that the unavailable differential opera- 
tor D xy corresponding to the microsimulator (|15|) does 
indeed possess the invariance property (0) for constants 
p = 3.0 and a = —2.0. 

To determine the self-similar shape, the template con- 
dition for the x direction was chosen as uj{— 2.266, oo, t) = 
0.4, i.e., the x-coordinates corresponding to the renormal- 
ized marginal CDF lujj = 0.4 always has the same value, 
—2.266cm. For the analytical solution (|22fl and our tem- 
plate c = 0.3125sec _1 and the corresponding std.'s are 
ax = 4\/5cm and cry = 12.8vT5/3cm, respectively. 

A uniform particle CDF over the space domain 
(—10cm, 10cm) x (—10cm, 10cm) is used as the initial 
condition for our equation-free fixed point algorithm. Di- 
rect iteration is used to solve the fixed point of equation 
1(11)1. The time interval V is 150Ai. The number of 
conditional ICDFs is 20 [M = 20) and the basis for the 
both types of ICDF is the shifted Legendre polynomials 
of order up to and including 5 (P = 5). In this simula- 
tion, 200 copies of 5000 particle positions are generated 
according to the ICDF mode coefficients at the beginning 
of each iteration and let to evolve for time T". The mode 
coefficients at the end of each iteration are obtained by 
averaging over these 200 copies. After the 2nd, 4th and 
6th iterations, mode coefficients of rcnormalizcd ICDFs 
are used to generate particle positions, out of which the 
CDFs are computed and plotted respectively in Fig. [3] 
It can be seen that the renormalized solutions reach a 
steady state quickly. 

To verify the self-similar shape of the solution, the 
std.'s and correlation of the simulated self-similar shape 
are compared with those of the analytical self-similar 
shape. The std.'s and correlations of the rescaled CDFs 
are calculated via the ensemble particle positions corre- 
sponding to these CDFs. The comparison is shown in 
Fig. 0] The std. 's and correlations of rescaled CDFs ap- 



t(sec) A(t) A t (t) 

~ 1.00000 

1.5 1.20872 0.13915 

3.0 1.38243 0.11581 



TABLE II: Values of the rescaling variable A(t) in Eqn.ltTSt 

proach those of the analytical self-similar shape, which 
means that the rescaled CDF quickly approaches a mem- 
ber of the family of theoretical self-similar shapes ex- 
pressed by Equation H22|) . 

As the renormalized CDF uj(u, v, t) reaches the steady 
state, we can set this CDF as the initial condition and 
run the microscale dynamics (|16f) for two more loops 
with t\ = 150At and ti = 300Ai. The rescaling vari- 
able A(t) is listed in Table HU Note that A(t) = 1 at 
t = 0. By Equation (|13(l . the similarity exponent a is 
approximately estimated as 0.461, which is within 8% of 
the theoretical value 1/2. 

We now study the effect of varying templates and re- 
porting horizons in the fixed-point operator <E>t', on our 
results; different templates should give different members 
in the family of Eqn. JSJ, while T 1 should not affect the 
computed fixed points. 

Four combinations of template condition and evolution 
times including the one above are investigated: 

1. w(-2.266,oo,i) = 0.4, T' = 150Ai; 

2. w(-2.266, oo, i) = 0.4, T 1 = 250Ai; 

3. w(-0.227,co,i) = 0.4, T' = 150Ai; (23) 

4. w(-0.227,oo,i) = 0.4, T' = 250A£. 

The iterative values of std.'s and correlation for the 
four cases are shown in figures 0] and [S] Comparison 
with theoretical calculations shows that variation of tem- 
plates and evolution times indeed does not cause devia- 
tion of the converged rescaled CDF from the family of 
self-similar solutions. 



VI. CONCLUSIONS 

We presented a coarse dynamic renormalization tech- 
nique, based on a coarse time-stepper using marginal 
and conditional ICDFs as coarse-scale observables, for 
the computer-assisted analysis of multidimensional self- 
similar particle systems. For coarse-grained differential 
operators that possess a scale-invariance property of the 
type arising in our numerical example, a single template 
condition is necessary for rescaling the ICDFs. In our ex- 
ample, We confirmed that the steady-state shape of the 
renormalized CDFs of particle positions is not affected 
by the template condition and evolution time interval in 
our coarse time-stepper. 

The equation-free dynamic renormalization technique 
was applied only to a simple 2D Brownian particle sys- 
tem in this paper. In our computations we had already 
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FIG. 4: Computational vs. theoretical values for the standard 
deviations and the correlation of the self-similar shapes. Case 
1, 2 in (23) (see text). 



FIG. 3: Renormalized CDFs ; top left: initial CDF, top right: 
renormalized CDF after 2nd iteration, bottom left: renor- 
malized CDF after 4th iteration, bottom right: renormalized 
CDF after 6th iteration (see text). 



factored out translational invariance; we also knew the 
particular coordinate system (the x and y axis) in which 
the macroscopic differential operator is scale invariant. 
Finding this information, if unknown, can be incorpo- 
rated in the test procedure for scale invariance. For 
higher-dimensional diffusive particle systems, this tech- 
nique may still be used. In these cases, conditional 
fCDFs in third or higher dimensions have to be utilized, 
thus rendering the computations slightly more compli- 
cated, yet still qualitatively similar to the ones presented 
here. Using appropriately selected template conditions, 
this technique can also be used to approximate coarse- 
grained asymptotically self-similar shapes; this is the case 
for our Couette example when the particles diffuse in both 
the x and y directions (work in progress). Several other 
examples of coarse self-similarity (e.g. models of glassy 
dynamics, KPZ-type evolution of interfaces, dynamics of 
energy spectra in randomly forced PDEs) are also being 
explored. 
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computational, case 3 in (23) 
computational, case 4 in (23) 
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FIG. 5: Computational vs. theoretical values for the standard 
deviations and the correlation of the self-similar shapes. Case 
3, 4 in (23) (see text). 
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